
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!
 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Linear Algebra Data and Routines File
! 
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



MODULE aqchem_LinearAlgebra

  USE aqchem_Parameters
  USE aqchem_JacobianSP

  IMPLICIT NONE

CONTAINS


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! SPARSE_UTIL - SPARSE utility functions
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecomp( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER  :: IER
      REAL(kind=dp) :: JVS(LU_NONZERO), W(NVAR), a
      INTEGER  :: k, kk, j, jj

      a = 0. ! mz_rs_20050606
      IER = 0
      DO k=1,NVAR
        ! mz_rs_20050606: don't check if real value == 0
        ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecomp


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplx( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization, complex
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: IER
      DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a
      REAL(kind=dp)  :: b = 0.0
      INTEGER        :: k, kk, j, jj

      IER = 0
      DO k=1,NVAR
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecompCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Sparse LU factorization, complex
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       :: IER
      REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO) 
      REAL(kind=dp) :: WR(NVAR), WI(NVAR), ar, ai, den
      INTEGER       :: k, kk, j, jj

      IER = 0
      ar  = 0.0
      DO k=1,NVAR
        IF (  ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. &
              ( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) )  THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              WR( LU_ICOL(kk) ) = JVSR(kk)
              WI( LU_ICOL(kk) ) = JVSI(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            den = JVSR(LU_DIAG(j))**2 + JVSI(LU_DIAG(j))**2
            ar = -(WR(j)*JVSR(LU_DIAG(j)) + WI(j)*JVSI(LU_DIAG(j)))/den
            ai = -(WI(j)*JVSR(LU_DIAG(j)) - WR(j)*JVSI(LU_DIAG(j)))/den
            WR(j) = -ar
            WI(j) = -ai
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + ar*JVSR(jj) - ai*JVSI(jj)
               WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + ar*JVSI(jj) + ai*JVSR(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVSR(kk) = WR( LU_ICOL(kk) )
            JVSI(kk) = WI( LU_ICOL(kk) )
         END DO
      END DO

END SUBROUTINE KppDecompCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER  :: i, j
      REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       :: i, j
      REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveCmplx

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       ::  i, j
      REAL(kind=dp) ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             XR(i) = XR(i) - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
             XI(i) = XI(i) - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
         END DO  
      END DO

      DO i=NVAR,1,-1
        sumr = XR(i); sumi = XI(i)
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
            sumr = sumr - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
            sumi = sumi - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
        END DO
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (sumr*JVSR(LU_DIAG(i)) + sumi*JVSI(LU_DIAG(i)))/den
        XI(i) = (sumi*JVSR(LU_DIAG(i)) - sumr*JVSI(LU_DIAG(i)))/den
      END DO
      
END SUBROUTINE KppSolveCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve transpose subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       ::  i, j
      REAL(kind=dp) ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den

      DO i=1,NVAR
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (XR(i)*JVSR(LU_DIAG(i)) + XI(i)*JVSI(LU_DIAG(i)))/den
        XI(i) = (XI(i)*JVSR(LU_DIAG(i)) - XR(i)*JVSI(LU_DIAG(i)))/den
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplxR


!
! Next few commented subroutines perform sparse big linear algebra
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppDecompBig( JVS, IP, IER )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse LU factorization
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER  :: IP3(3), IER, IP(3,NVAR)
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3)
!      INTEGER  :: k, kk, j, jj
!
!      a = 0.0d0
!      IER = 0
!      DO k=1,NVAR
!        DO kk = LU_CROW(k), LU_CROW(k+1)-1
!              W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk)
!        END DO
!        DO kk = LU_CROW(k), LU_DIAG(k)-1
!            j = LU_ICOL(kk)
!            E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) )
!            ! CALL DGETRF(3,3,E,3,IP3,IER) 
!            CALL FAC3(E,IP3,IER)
!            IF ( IER /= 0 )  RETURN
!            ! a = W(j) / JVS( LU_DIAG(j) )
!            a(1:3,1:3) = W( 1:3,1:3,j )
!            ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER) 
!            CALL SOL3('N',E,IP3,a(1,1))
!            CALL SOL3('N',E,IP3,a(1,2))
!            CALL SOL3('N',E,IP3,a(1,3))
!            W(1:3,1:3,j) = a(1:3,1:3)
!            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
!               W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) &
!                        - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) )
!            END DO
!         END DO
!         DO kk = LU_CROW(k), LU_CROW(k+1)-1
!            JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) )
!         END DO
!      END DO
!
!      DO k=1,NVAR
!         ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER)
!         ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER)
!         CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER)
!         IF ( IER /= 0 )  RETURN
!      END DO 
!      
!END SUBROUTINE KppDecompBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBig( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse solve subroutine using indirect addressing
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER  :: i, j, k, m, IP3(3), IP(3,NVAR), IER
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3)
!
!      DO i=1,NVAR
!        DO j = LU_CROW(i), LU_DIAG(i)-1 
!          !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!               X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO  
!      END DO
!
!      DO i=NVAR,1,-1
!        sum(1:3) = X(1:3,i);
!        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
!          !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!               sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO
!        ! X(i) = sum/JVS(LU_DIAG(i));
!        ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0) 
!        ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum)
!        CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum)
!        X(1:3,i) = sum(1:3)
!      END DO
!      
!END SUBROUTINE KppSolveBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBigTR( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Big sparse transpose solve using indirect addressing
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER       :: i, j, k, m, IP(3,NVAR)
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR)
!
!      DO i=1,NVAR
!        ! X(i) = X(i)/JVS(LU_DIAG(i))
!        CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i))
!        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
!          !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !    - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!               X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!        END DO
!      END DO
!
!      DO i=NVAR, 1, -1
!        DO j=LU_CROW(i),LU_DIAG(i)-1
!          !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !   - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!               X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!        END DO
!      END DO
!      
!END SUBROUTINE KppSolveBigTR
!
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE FAC3(A,IPVT,INFO)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     FAC3 FACTORS THE MATRIX A (3,3) BY
!!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
!!     LINPACK - LIKE 
!!
!!     Remove comments to perform pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!
!      REAL(kind=dp) :: A(3,3)
!      INTEGER       :: IPVT(3),INFO
!!      INTEGER       :: L
!!      REAL(kind=dp) :: t, dmax, da, TMP(3)
!      REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0
!
!      info = 0
!!      t = TINY(da)
!!      
!!      da = ABS(A(1,1)); L = 1
!!      IF ( ABS(A(2,1))>da ) THEN
!!        da = ABS(A(2,1)); L = 2
!!        IF ( ABS(A(3,1))>da ) THEN
!!          L = 3
!!        END IF  
!!      END IF  
!!      IPVT(1)  = L
!!      IF (L /=1 ) THEN
!!         TMP(1:3) = A(L,1:3)
!!         A(L,1:3) = A(1,1:3)
!!         A(1,1:3) = TMP(1:3)
!!      END IF
!!      IF (ABS(A(1,1)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!
!      A(2,1) = A(2,1)/A(1,1)
!      A(2,2) = A(2,2) - A(2,1)*A(1,2)
!      A(2,3) = A(2,3) - A(2,1)*A(1,3)
!      A(3,1) = A(3,1)/A(1,1)
!      A(3,2) = A(3,2) - A(3,1)*A(1,2)
!      A(3,3) = A(3,3) - A(3,1)*A(1,3)
!      
!!      IPVT(2)  = 2
!!      IF (ABS(A(3,2))>ABS(A(2,2))) THEN
!!         IPVT(2)  = 3
!!         TMP(2:3) = A(3,2:3)
!!         A(3,2:3) = A(2,2:3)
!!         A(2,2:3) = TMP(2:3)
!!      END IF
!!      IF (ABS(A(2,2)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!      
!      A(3,2)   = A(3,2)/A(2,2)
!      A(3,3)   = A(3,3) - A(3,2)*A(2,3)
!      IPVT(3)  = 3
!      
!END SUBROUTINE FAC3
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE SOL3(Trans,A,IPVT,b)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     SOL3 solves the system 3x3
!!     A * x = b  or  trans(a) * x = b
!!     using the factors computed by WGEFA.
!!
!!     Trans      = 'N'   to solve  A*x = b ,
!!                = 'T'   to solve  transpose(A)*x = b
!!     LINPACK - LIKE 
!!
!!     Remove comments to use pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!      CHARACTER     :: Trans
!      REAL(kind=dp) :: a(3,3),b(3)
!      INTEGER       :: IPVT(3)
!!      INTEGER       :: L
!!      REAL(kind=dp) :: TMP
!      
!      SELECT CASE (Trans)
!
!      CASE ('n','N')  !  Solve  A * x = b
!
!!     Solve  L*y = b
!!         L = IPVT(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!         b(2) = b(2)-A(2,1)*b(1)
!         b(3) = b(3)-A(3,1)*b(1)
!         
!!         L = IPVT(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(3) = b(3)-A(3,2)*b(2)
!
!!     Solve  U*x = y
!         b(3) = b(3)/A(3,3)
!         b(2) = (b(2)-A(2,3)*b(3))/A(2,2)
!         b(1) = (b(1)-A(1,3)*b(3)-A(1,2)*b(2))/A(1,1)
!      
!      
!      CASE ('t','T')  !  Solve transpose(A) * x = b
!
!!      Solve transpose(U)*y = b
!         b(1) = b(1)/A(1,1)
!         b(2) = (b(2)-A(1,2)*b(1))/A(2,2)
!         b(3) = (b(3)-A(1,3)*b(1)-A(2,3)*b(2))/A(3,3)
!
!!      Solve transpose(L)*x = y
!         b(2) = b(2)-A(3,2)*b(3)
!!         L = ipvt(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(1) = b(1)-A(3,1)*b(3)-A(2,1)*b(2)
!!         L = ipvt(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!   
!      END SELECT
!
!END SUBROUTINE SOL3

! End of SPARSE_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! KppSolve - sparse back substitution
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      X         - Vector for variables
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE KppSolve ( JVS, X )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! X - Vector for variables
  REAL(kind=dp) :: X(NVAR)

  X(9) = X(9)-JVS(9)*X(8)
  X(10) = X(10)-JVS(11)*X(6)
  X(12) = X(12)-JVS(14)*X(7)
  X(18) = X(18)-JVS(21)*X(4)
  X(37) = X(37)-JVS(74)*X(9)
  X(38) = X(38)-JVS(76)*X(13)
  X(39) = X(39)-JVS(78)*X(14)
  X(40) = X(40)-JVS(80)*X(10)
  X(41) = X(41)-JVS(82)*X(11)
  X(42) = X(42)-JVS(84)*X(12)
  X(44) = X(44)-JVS(88)*X(15)
  X(45) = X(45)-JVS(90)*X(16)
  X(46) = X(46)-JVS(92)*X(17)
  X(47) = X(47)-JVS(94)*X(18)
  X(76) = X(76)-JVS(173)*X(75)
  X(78) = X(78)-JVS(178)*X(77)
  X(105) = X(105)-JVS(237)*X(79)
  X(106) = X(106)-JVS(240)*X(94)-JVS(241)*X(105)
  X(107) = X(107)-JVS(244)*X(86)
  X(108) = X(108)-JVS(248)*X(88)
  X(109) = X(109)-JVS(252)*X(2)-JVS(253)*X(108)
  X(111) = X(111)-JVS(259)*X(93)
  X(112) = X(112)-JVS(263)*X(3)-JVS(264)*X(111)
  X(113) = X(113)-JVS(267)*X(95)
  X(115) = X(115)-JVS(273)*X(114)
  X(116) = X(116)-JVS(277)*X(101)-JVS(278)*X(102)-JVS(279)*X(103)
  X(119) = X(119)-JVS(294)*X(84)-JVS(295)*X(117)-JVS(296)*X(118)
  X(120) = X(120)-JVS(302)*X(99)
  X(121) = X(121)-JVS(305)*X(91)-JVS(306)*X(120)
  X(122) = X(122)-JVS(310)*X(92)
  X(124) = X(124)-JVS(318)*X(110)-JVS(319)*X(123)
  X(125) = X(125)-JVS(325)*X(87)-JVS(326)*X(115)
  X(126) = X(126)-JVS(337)*X(74)
  X(127) = X(127)-JVS(342)*X(98)
  X(132) = X(132)-JVS(365)*X(129)-JVS(366)*X(131)
  X(133) = X(133)-JVS(372)*X(113)
  X(134) = X(134)-JVS(377)*X(83)-JVS(378)*X(133)
  X(135) = X(135)-JVS(382)*X(80)-JVS(383)*X(134)
  X(136) = X(136)-JVS(388)*X(133)-JVS(389)*X(134)-JVS(390)*X(135)
  X(137) = X(137)-JVS(394)*X(106)-JVS(395)*X(135)-JVS(396)*X(136)
  X(139) = X(139)-JVS(405)*X(105)
  X(140) = X(140)-JVS(410)*X(139)
  X(141) = X(141)-JVS(414)*X(113)-JVS(415)*X(136)-JVS(416)*X(137)-JVS(417)*X(139)-JVS(418)*X(140)
  X(142) = X(142)-JVS(422)*X(90)-JVS(423)*X(138)
  X(143) = X(143)-JVS(429)*X(104)-JVS(430)*X(120)-JVS(431)*X(127)
  X(144) = X(144)-JVS(437)*X(81)-JVS(438)*X(120)
  X(145) = X(145)-JVS(444)*X(13)-JVS(445)*X(14)-JVS(446)*X(100)
  X(146) = X(146)-JVS(452)*X(85)
  X(147) = X(147)-JVS(459)*X(1)-JVS(460)*X(107)-JVS(461)*X(116)-JVS(462)*X(127)
  X(148) = X(148)-JVS(470)*X(144)
  X(149) = X(149)-JVS(477)*X(135)-JVS(478)*X(136)-JVS(479)*X(144)-JVS(480)*X(148)
  X(150) = X(150)-JVS(487)*X(97)-JVS(488)*X(130)
  X(152) = X(152)-JVS(501)*X(126)-JVS(502)*X(130)-JVS(503)*X(145)-JVS(504)*X(150)-JVS(505)*X(151)
  X(153) = X(153)-JVS(514)*X(108)-JVS(515)*X(109)-JVS(516)*X(131)-JVS(517)*X(132)-JVS(518)*X(137)-JVS(519)*X(138)&
             &-JVS(520)*X(141)-JVS(521)*X(142)-JVS(522)*X(144)-JVS(523)*X(147)-JVS(524)*X(148)-JVS(525)*X(149)-JVS(526)&
             &*X(150)-JVS(527)*X(151)-JVS(528)*X(152)
  X(154) = X(154)-JVS(538)*X(5)-JVS(539)*X(13)-JVS(540)*X(14)-JVS(541)*X(76)-JVS(542)*X(78)-JVS(543)*X(116)-JVS(544)&
             &*X(121)-JVS(545)*X(122)-JVS(546)*X(124)-JVS(547)*X(127)-JVS(548)*X(128)-JVS(549)*X(143)-JVS(550)*X(145)&
             &-JVS(551)*X(146)-JVS(552)*X(147)-JVS(553)*X(150)-JVS(554)*X(151)-JVS(555)*X(152)
  X(155) = X(155)-JVS(564)*X(82)-JVS(565)*X(105)-JVS(566)*X(106)-JVS(567)*X(113)-JVS(568)*X(120)-JVS(569)*X(133)&
             &-JVS(570)*X(134)-JVS(571)*X(137)-JVS(572)*X(138)-JVS(573)*X(139)-JVS(574)*X(140)-JVS(575)*X(141)-JVS(576)&
             &*X(142)-JVS(577)*X(143)-JVS(578)*X(146)-JVS(579)*X(148)-JVS(580)*X(149)-JVS(581)*X(150)-JVS(582)*X(151)&
             &-JVS(583)*X(153)-JVS(584)*X(154)
  X(156) = X(156)-JVS(592)*X(13)-JVS(593)*X(14)-JVS(594)*X(107)-JVS(595)*X(110)-JVS(596)*X(111)-JVS(597)*X(112)-JVS(598)&
             &*X(114)-JVS(599)*X(115)-JVS(600)*X(121)-JVS(601)*X(122)-JVS(602)*X(124)-JVS(603)*X(125)-JVS(604)*X(126)&
             &-JVS(605)*X(127)-JVS(606)*X(128)-JVS(607)*X(129)-JVS(608)*X(130)-JVS(609)*X(131)-JVS(610)*X(132)-JVS(611)&
             &*X(133)-JVS(612)*X(134)-JVS(613)*X(135)-JVS(614)*X(136)-JVS(615)*X(137)-JVS(616)*X(138)-JVS(617)*X(139)&
             &-JVS(618)*X(140)-JVS(619)*X(141)-JVS(620)*X(142)-JVS(621)*X(143)-JVS(622)*X(145)-JVS(623)*X(146)-JVS(624)&
             &*X(147)-JVS(625)*X(148)-JVS(626)*X(149)-JVS(627)*X(150)-JVS(628)*X(151)-JVS(629)*X(152)-JVS(630)*X(153)&
             &-JVS(631)*X(154)-JVS(632)*X(155)
  X(157) = X(157)-JVS(639)*X(13)-JVS(640)*X(14)-JVS(641)*X(120)-JVS(642)*X(121)-JVS(643)*X(122)-JVS(644)*X(127)-JVS(645)&
             &*X(128)-JVS(646)*X(143)-JVS(647)*X(144)-JVS(648)*X(145)-JVS(649)*X(146)-JVS(650)*X(148)-JVS(651)*X(149)&
             &-JVS(652)*X(150)-JVS(653)*X(151)-JVS(654)*X(153)-JVS(655)*X(154)-JVS(656)*X(155)-JVS(657)*X(156)
  X(158) = X(158)-JVS(663)*X(89)-JVS(664)*X(145)-JVS(665)*X(151)-JVS(666)*X(152)-JVS(667)*X(155)-JVS(668)*X(156)&
             &-JVS(669)*X(157)
  X(159) = X(159)-JVS(674)*X(96)-JVS(675)*X(126)-JVS(676)*X(145)-JVS(677)*X(147)-JVS(678)*X(150)-JVS(679)*X(151)&
             &-JVS(680)*X(152)-JVS(681)*X(154)-JVS(682)*X(155)-JVS(683)*X(156)-JVS(684)*X(157)-JVS(685)*X(158)
  X(160) = X(160)-JVS(689)*X(13)-JVS(690)*X(14)-JVS(691)*X(144)-JVS(692)*X(148)-JVS(693)*X(149)-JVS(694)*X(153)-JVS(695)&
             &*X(154)-JVS(696)*X(155)-JVS(697)*X(156)-JVS(698)*X(157)-JVS(699)*X(158)-JVS(700)*X(159)
  X(161) = X(161)-JVS(703)*X(105)-JVS(704)*X(106)-JVS(705)*X(113)-JVS(706)*X(123)-JVS(707)*X(126)-JVS(708)*X(127)&
             &-JVS(709)*X(129)-JVS(710)*X(131)-JVS(711)*X(132)-JVS(712)*X(133)-JVS(713)*X(134)-JVS(714)*X(135)-JVS(715)&
             &*X(136)-JVS(716)*X(137)-JVS(717)*X(138)-JVS(718)*X(139)-JVS(719)*X(140)-JVS(720)*X(141)-JVS(721)*X(142)&
             &-JVS(722)*X(143)-JVS(723)*X(146)-JVS(724)*X(147)-JVS(725)*X(148)-JVS(726)*X(149)-JVS(727)*X(150)-JVS(728)&
             &*X(151)-JVS(729)*X(152)-JVS(730)*X(153)-JVS(731)*X(154)-JVS(732)*X(155)-JVS(733)*X(156)-JVS(734)*X(157)&
             &-JVS(735)*X(158)-JVS(736)*X(159)-JVS(737)*X(160)
  X(161) = X(161)/JVS(738)
  X(160) = (X(160)-JVS(702)*X(161))/(JVS(701))
  X(159) = (X(159)-JVS(687)*X(160)-JVS(688)*X(161))/(JVS(686))
  X(158) = (X(158)-JVS(671)*X(159)-JVS(672)*X(160)-JVS(673)*X(161))/(JVS(670))
  X(157) = (X(157)-JVS(659)*X(158)-JVS(660)*X(159)-JVS(661)*X(160)-JVS(662)*X(161))/(JVS(658))
  X(156) = (X(156)-JVS(634)*X(157)-JVS(635)*X(158)-JVS(636)*X(159)-JVS(637)*X(160)-JVS(638)*X(161))/(JVS(633))
  X(155) = (X(155)-JVS(586)*X(156)-JVS(587)*X(157)-JVS(588)*X(158)-JVS(589)*X(159)-JVS(590)*X(160)-JVS(591)*X(161))&
             &/(JVS(585))
  X(154) = (X(154)-JVS(557)*X(155)-JVS(558)*X(156)-JVS(559)*X(157)-JVS(560)*X(158)-JVS(561)*X(159)-JVS(562)*X(160)&
             &-JVS(563)*X(161))/(JVS(556))
  X(153) = (X(153)-JVS(530)*X(154)-JVS(531)*X(155)-JVS(532)*X(156)-JVS(533)*X(157)-JVS(534)*X(158)-JVS(535)*X(159)&
             &-JVS(536)*X(160)-JVS(537)*X(161))/(JVS(529))
  X(152) = (X(152)-JVS(507)*X(155)-JVS(508)*X(156)-JVS(509)*X(157)-JVS(510)*X(158)-JVS(511)*X(159)-JVS(512)*X(160)&
             &-JVS(513)*X(161))/(JVS(506))
  X(151) = (X(151)-JVS(496)*X(155)-JVS(497)*X(156)-JVS(498)*X(158)-JVS(499)*X(159)-JVS(500)*X(161))/(JVS(495))
  X(150) = (X(150)-JVS(490)*X(151)-JVS(491)*X(155)-JVS(492)*X(156)-JVS(493)*X(157)-JVS(494)*X(159))/(JVS(489))
  X(149) = (X(149)-JVS(482)*X(153)-JVS(483)*X(156)-JVS(484)*X(157)-JVS(485)*X(160)-JVS(486)*X(161))/(JVS(481))
  X(148) = (X(148)-JVS(472)*X(149)-JVS(473)*X(153)-JVS(474)*X(157)-JVS(475)*X(160)-JVS(476)*X(161))/(JVS(471))
  X(147) = (X(147)-JVS(464)*X(150)-JVS(465)*X(152)-JVS(466)*X(154)-JVS(467)*X(156)-JVS(468)*X(157)-JVS(469)*X(158))&
             &/(JVS(463))
  X(146) = (X(146)-JVS(454)*X(151)-JVS(455)*X(155)-JVS(456)*X(156)-JVS(457)*X(157)-JVS(458)*X(161))/(JVS(453))
  X(145) = (X(145)-JVS(448)*X(156)-JVS(449)*X(157)-JVS(450)*X(158)-JVS(451)*X(159))/(JVS(447))
  X(144) = (X(144)-JVS(440)*X(148)-JVS(441)*X(149)-JVS(442)*X(157)-JVS(443)*X(160))/(JVS(439))
  X(143) = (X(143)-JVS(433)*X(154)-JVS(434)*X(155)-JVS(435)*X(157)-JVS(436)*X(161))/(JVS(432))
  X(142) = (X(142)-JVS(425)*X(148)-JVS(426)*X(149)-JVS(427)*X(156)-JVS(428)*X(161))/(JVS(424))
  X(141) = (X(141)-JVS(420)*X(156)-JVS(421)*X(161))/(JVS(419))
  X(140) = (X(140)-JVS(412)*X(156)-JVS(413)*X(161))/(JVS(411))
  X(139) = (X(139)-JVS(407)*X(140)-JVS(408)*X(156)-JVS(409)*X(161))/(JVS(406))
  X(138) = (X(138)-JVS(402)*X(142)-JVS(403)*X(156)-JVS(404)*X(161))/(JVS(401))
  X(137) = (X(137)-JVS(398)*X(141)-JVS(399)*X(156)-JVS(400)*X(161))/(JVS(397))
  X(136) = (X(136)-JVS(392)*X(156)-JVS(393)*X(161))/(JVS(391))
  X(135) = (X(135)-JVS(385)*X(136)-JVS(386)*X(156)-JVS(387)*X(161))/(JVS(384))
  X(134) = (X(134)-JVS(380)*X(156)-JVS(381)*X(161))/(JVS(379))
  X(133) = (X(133)-JVS(374)*X(134)-JVS(375)*X(156)-JVS(376)*X(161))/(JVS(373))
  X(132) = (X(132)-JVS(368)*X(137)-JVS(369)*X(141)-JVS(370)*X(156)-JVS(371)*X(161))/(JVS(367))
  X(131) = (X(131)-JVS(362)*X(132)-JVS(363)*X(156)-JVS(364)*X(161))/(JVS(361))
  X(130) = (X(130)-JVS(357)*X(150)-JVS(358)*X(151)-JVS(359)*X(156)-JVS(360)*X(159))/(JVS(356))
  X(129) = (X(129)-JVS(352)*X(132)-JVS(353)*X(137)-JVS(354)*X(156)-JVS(355)*X(161))/(JVS(351))
  X(128) = (X(128)-JVS(347)*X(143)-JVS(348)*X(155)-JVS(349)*X(156)-JVS(350)*X(157))/(JVS(346))
  X(127) = (X(127)-JVS(344)*X(154)-JVS(345)*X(157))/(JVS(343))
  X(126) = (X(126)-JVS(339)*X(152)-JVS(340)*X(156)-JVS(341)*X(161))/(JVS(338))
  X(125) = (X(125)-JVS(328)*X(129)-JVS(329)*X(131)-JVS(330)*X(132)-JVS(331)*X(133)-JVS(332)*X(134)-JVS(333)*X(138)&
             &-JVS(334)*X(142)-JVS(335)*X(156)-JVS(336)*X(161))/(JVS(327))
  X(124) = (X(124)-JVS(321)*X(127)-JVS(322)*X(143)-JVS(323)*X(154)-JVS(324)*X(156))/(JVS(320))
  X(123) = (X(123)-JVS(315)*X(127)-JVS(316)*X(143)-JVS(317)*X(154))/(JVS(314))
  X(122) = (X(122)-JVS(312)*X(156)-JVS(313)*X(157))/(JVS(311))
  X(121) = (X(121)-JVS(308)*X(156)-JVS(309)*X(157))/(JVS(307))
  X(120) = (X(120)-JVS(304)*X(157))/(JVS(303))
  X(119) = (X(119)-JVS(298)*X(153)-JVS(299)*X(157)-JVS(300)*X(160)-JVS(301)*X(161))/(JVS(297))
  X(118) = (X(118)-JVS(289)*X(119)-JVS(290)*X(153)-JVS(291)*X(157)-JVS(292)*X(160)-JVS(293)*X(161))/(JVS(288))
  X(117) = (X(117)-JVS(284)*X(118)-JVS(285)*X(119)-JVS(286)*X(153)-JVS(287)*X(161))/(JVS(283))
  X(116) = (X(116)-JVS(281)*X(147)-JVS(282)*X(154))/(JVS(280))
  X(115) = (X(115)-JVS(275)*X(125)-JVS(276)*X(156))/(JVS(274))
  X(114) = (X(114)-JVS(271)*X(115)-JVS(272)*X(156))/(JVS(270))
  X(113) = (X(113)-JVS(269)*X(161))/(JVS(268))
  X(112) = (X(112)-JVS(266)*X(156))/(JVS(265))
  X(111) = (X(111)-JVS(261)*X(112)-JVS(262)*X(156))/(JVS(260))
  X(110) = (X(110)-JVS(257)*X(124)-JVS(258)*X(156))/(JVS(256))
  X(109) = (X(109)-JVS(255)*X(153))/(JVS(254))
  X(108) = (X(108)-JVS(250)*X(109)-JVS(251)*X(153))/(JVS(249))
  X(107) = (X(107)-JVS(246)*X(147)-JVS(247)*X(156))/(JVS(245))
  X(106) = (X(106)-JVS(243)*X(161))/(JVS(242))
  X(105) = (X(105)-JVS(239)*X(161))/(JVS(238))
  X(104) = (X(104)-JVS(233)*X(120)-JVS(234)*X(127)-JVS(235)*X(157)-JVS(236)*X(161))/(JVS(232))
  X(103) = (X(103)-JVS(230)*X(116)-JVS(231)*X(154))/(JVS(229))
  X(102) = (X(102)-JVS(228)*X(116))/(JVS(227))
  X(101) = (X(101)-JVS(226)*X(116))/(JVS(225))
  X(100) = (X(100)-JVS(224)*X(145))/(JVS(223))
  X(99) = (X(99)-JVS(222)*X(120))/(JVS(221))
  X(98) = (X(98)-JVS(220)*X(127))/(JVS(219))
  X(97) = (X(97)-JVS(218)*X(150))/(JVS(217))
  X(96) = (X(96)-JVS(216)*X(159))/(JVS(215))
  X(95) = (X(95)-JVS(214)*X(113))/(JVS(213))
  X(94) = (X(94)-JVS(212)*X(106))/(JVS(211))
  X(93) = (X(93)-JVS(210)*X(111))/(JVS(209))
  X(92) = (X(92)-JVS(208)*X(122))/(JVS(207))
  X(91) = (X(91)-JVS(206)*X(121))/(JVS(205))
  X(90) = (X(90)-JVS(204)*X(142))/(JVS(203))
  X(89) = (X(89)-JVS(202)*X(158))/(JVS(201))
  X(88) = (X(88)-JVS(200)*X(108))/(JVS(199))
  X(87) = (X(87)-JVS(198)*X(125))/(JVS(197))
  X(86) = (X(86)-JVS(196)*X(107))/(JVS(195))
  X(85) = (X(85)-JVS(194)*X(146))/(JVS(193))
  X(84) = (X(84)-JVS(192)*X(119))/(JVS(191))
  X(83) = (X(83)-JVS(190)*X(134))/(JVS(189))
  X(82) = (X(82)-JVS(188)*X(155))/(JVS(187))
  X(81) = (X(81)-JVS(186)*X(144))/(JVS(185))
  X(80) = (X(80)-JVS(184)*X(135))/(JVS(183))
  X(79) = (X(79)-JVS(182)*X(105))/(JVS(181))
  X(78) = (X(78)-JVS(180)*X(154))/(JVS(179))
  X(77) = (X(77)-JVS(177)*X(78))/(JVS(176))
  X(76) = (X(76)-JVS(175)*X(154))/(JVS(174))
  X(75) = (X(75)-JVS(172)*X(76))/(JVS(171))
  X(74) = (X(74)-JVS(170)*X(126))/(JVS(169))
  X(73) = (X(73)-JVS(168)*X(120))/(JVS(167))
  X(72) = (X(72)-JVS(166)*X(127))/(JVS(165))
  X(71) = (X(71)-JVS(162)*X(117)-JVS(163)*X(118)-JVS(164)*X(119))/(JVS(161))
  X(70) = (X(70)-JVS(159)*X(139)-JVS(160)*X(140))/(JVS(158))
  X(69) = (X(69)-JVS(157)*X(105))/(JVS(156))
  X(68) = (X(68)-JVS(154)*X(151)-JVS(155)*X(155))/(JVS(153))
  X(67) = (X(67)-JVS(150)*X(144)-JVS(151)*X(148)-JVS(152)*X(149))/(JVS(149))
  X(66) = (X(66)-JVS(147)*X(135)-JVS(148)*X(136))/(JVS(146))
  X(65) = (X(65)-JVS(144)*X(133)-JVS(145)*X(134))/(JVS(143))
  X(64) = (X(64)-JVS(140)*X(129)-JVS(141)*X(131)-JVS(142)*X(132))/(JVS(139))
  X(63) = (X(63)-JVS(137)*X(137)-JVS(138)*X(141))/(JVS(136))
  X(62) = (X(62)-JVS(133)*X(76)-JVS(134)*X(78)-JVS(135)*X(154))/(JVS(132))
  X(61) = (X(61)-JVS(131)*X(62))/(JVS(130))
  X(60) = (X(60)-JVS(128)*X(76)-JVS(129)*X(78))/(JVS(127))
  X(59) = (X(59)-JVS(126)*X(60))/(JVS(125))
  X(58) = (X(58)-JVS(124)*X(78))/(JVS(123))
  X(57) = (X(57)-JVS(122)*X(76))/(JVS(121))
  X(56) = (X(56)-JVS(118)*X(102)-JVS(119)*X(103)-JVS(120)*X(116))/(JVS(117))
  X(55) = (X(55)-JVS(116)*X(56))/(JVS(115))
  X(54) = (X(54)-JVS(114)*X(103))/(JVS(113))
  X(53) = (X(53)-JVS(110)*X(116)-JVS(111)*X(147)-JVS(112)*X(154))/(JVS(109))
  X(52) = (X(52)-JVS(107)*X(53)-JVS(108)*X(102))/(JVS(106))
  X(51) = (X(51)-JVS(105)*X(116))/(JVS(104))
  X(50) = (X(50)-JVS(102)*X(130)-JVS(103)*X(150))/(JVS(101))
  X(49) = (X(49)-JVS(99)*X(126)-JVS(100)*X(152))/(JVS(98))
  X(48) = (X(48)-JVS(97)*X(159))/(JVS(96))
  X(47) = X(47)/JVS(95)
  X(46) = X(46)/JVS(93)
  X(45) = X(45)/JVS(91)
  X(44) = X(44)/JVS(89)
  X(43) = (X(43)-JVS(87)*X(156))/(JVS(86))
  X(42) = X(42)/JVS(85)
  X(41) = X(41)/JVS(83)
  X(40) = X(40)/JVS(81)
  X(39) = X(39)/JVS(79)
  X(38) = X(38)/JVS(77)
  X(37) = X(37)/JVS(75)
  X(36) = (X(36)-JVS(73)*X(112))/(JVS(72))
  X(35) = (X(35)-JVS(71)*X(109))/(JVS(70))
  X(34) = (X(34)-JVS(69)*X(147))/(JVS(68))
  X(33) = (X(33)-JVS(67)*X(161))/(JVS(66))
  X(32) = (X(32)-JVS(65)*X(113))/(JVS(64))
  X(31) = (X(31)-JVS(63)*X(106))/(JVS(62))
  X(30) = (X(30)-JVS(61)*X(111))/(JVS(60))
  X(29) = (X(29)-JVS(57)*X(110)-JVS(58)*X(124)-JVS(59)*X(154))/(JVS(56))
  X(28) = (X(28)-JVS(55)*X(122))/(JVS(54))
  X(27) = (X(27)-JVS(53)*X(121))/(JVS(52))
  X(26) = (X(26)-JVS(50)*X(138)-JVS(51)*X(142))/(JVS(49))
  X(25) = (X(25)-JVS(48)*X(158))/(JVS(47))
  X(24) = (X(24)-JVS(46)*X(146))/(JVS(45))
  X(23) = (X(23)-JVS(44)*X(108))/(JVS(43))
  X(22) = (X(22)-JVS(40)*X(131)-JVS(41)*X(132)-JVS(42)*X(161))/(JVS(39))
  X(21) = (X(21)-JVS(35)*X(22)-JVS(36)*X(114)-JVS(37)*X(115)-JVS(38)*X(125))/(JVS(34))
  X(20) = (X(20)-JVS(33)*X(107))/(JVS(32))
  X(19) = (X(19)-JVS(24)*X(104)-JVS(25)*X(123)-JVS(26)*X(128)-JVS(27)*X(143)-JVS(28)*X(145)-JVS(29)*X(148)-JVS(30)&
            &*X(157)-JVS(31)*X(160))/(JVS(23))
  X(18) = X(18)/JVS(22)
  X(17) = X(17)/JVS(20)
  X(16) = X(16)/JVS(19)
  X(15) = X(15)/JVS(18)
  X(14) = X(14)/JVS(17)
  X(13) = X(13)/JVS(16)
  X(12) = X(12)/JVS(15)
  X(11) = X(11)/JVS(13)
  X(10) = X(10)/JVS(12)
  X(9) = X(9)/JVS(10)
  X(8) = X(8)/JVS(8)
  X(7) = X(7)/JVS(7)
  X(6) = X(6)/JVS(6)
  X(5) = X(5)/JVS(5)
  X(4) = X(4)/JVS(4)
  X(3) = X(3)/JVS(3)
  X(2) = X(2)/JVS(2)
  X(1) = X(1)/JVS(1)
      
END SUBROUTINE KppSolve

! End of KppSolve function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! KppSolveTR - sparse, transposed back substitution
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      X         - Vector for variables
!      XX        - Vector for output variables
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE KppSolveTR ( JVS, X, XX )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! X - Vector for variables
  REAL(kind=dp) :: X(NVAR)
! XX - Vector for output variables
  REAL(kind=dp) :: XX(NVAR)

  XX(1) = X(1)/JVS(1)
  XX(2) = X(2)/JVS(2)
  XX(3) = X(3)/JVS(3)
  XX(4) = X(4)/JVS(4)
  XX(5) = X(5)/JVS(5)
  XX(6) = X(6)/JVS(6)
  XX(7) = X(7)/JVS(7)
  XX(8) = X(8)/JVS(8)
  XX(9) = X(9)/JVS(10)
  XX(10) = X(10)/JVS(12)
  XX(11) = X(11)/JVS(13)
  XX(12) = X(12)/JVS(15)
  XX(13) = X(13)/JVS(16)
  XX(14) = X(14)/JVS(17)
  XX(15) = X(15)/JVS(18)
  XX(16) = X(16)/JVS(19)
  XX(17) = X(17)/JVS(20)
  XX(18) = X(18)/JVS(22)
  XX(19) = X(19)/JVS(23)
  XX(20) = X(20)/JVS(32)
  XX(21) = X(21)/JVS(34)
  XX(22) = (X(22)-JVS(35)*XX(21))/(JVS(39))
  XX(23) = X(23)/JVS(43)
  XX(24) = X(24)/JVS(45)
  XX(25) = X(25)/JVS(47)
  XX(26) = X(26)/JVS(49)
  XX(27) = X(27)/JVS(52)
  XX(28) = X(28)/JVS(54)
  XX(29) = X(29)/JVS(56)
  XX(30) = X(30)/JVS(60)
  XX(31) = X(31)/JVS(62)
  XX(32) = X(32)/JVS(64)
  XX(33) = X(33)/JVS(66)
  XX(34) = X(34)/JVS(68)
  XX(35) = X(35)/JVS(70)
  XX(36) = X(36)/JVS(72)
  XX(37) = X(37)/JVS(75)
  XX(38) = X(38)/JVS(77)
  XX(39) = X(39)/JVS(79)
  XX(40) = X(40)/JVS(81)
  XX(41) = X(41)/JVS(83)
  XX(42) = X(42)/JVS(85)
  XX(43) = X(43)/JVS(86)
  XX(44) = X(44)/JVS(89)
  XX(45) = X(45)/JVS(91)
  XX(46) = X(46)/JVS(93)
  XX(47) = X(47)/JVS(95)
  XX(48) = X(48)/JVS(96)
  XX(49) = X(49)/JVS(98)
  XX(50) = X(50)/JVS(101)
  XX(51) = X(51)/JVS(104)
  XX(52) = X(52)/JVS(106)
  XX(53) = (X(53)-JVS(107)*XX(52))/(JVS(109))
  XX(54) = X(54)/JVS(113)
  XX(55) = X(55)/JVS(115)
  XX(56) = (X(56)-JVS(116)*XX(55))/(JVS(117))
  XX(57) = X(57)/JVS(121)
  XX(58) = X(58)/JVS(123)
  XX(59) = X(59)/JVS(125)
  XX(60) = (X(60)-JVS(126)*XX(59))/(JVS(127))
  XX(61) = X(61)/JVS(130)
  XX(62) = (X(62)-JVS(131)*XX(61))/(JVS(132))
  XX(63) = X(63)/JVS(136)
  XX(64) = X(64)/JVS(139)
  XX(65) = X(65)/JVS(143)
  XX(66) = X(66)/JVS(146)
  XX(67) = X(67)/JVS(149)
  XX(68) = X(68)/JVS(153)
  XX(69) = X(69)/JVS(156)
  XX(70) = X(70)/JVS(158)
  XX(71) = X(71)/JVS(161)
  XX(72) = X(72)/JVS(165)
  XX(73) = X(73)/JVS(167)
  XX(74) = X(74)/JVS(169)
  XX(75) = X(75)/JVS(171)
  XX(76) = (X(76)-JVS(122)*XX(57)-JVS(128)*XX(60)-JVS(133)*XX(62)-JVS(172)*XX(75))/(JVS(174))
  XX(77) = X(77)/JVS(176)
  XX(78) = (X(78)-JVS(124)*XX(58)-JVS(129)*XX(60)-JVS(134)*XX(62)-JVS(177)*XX(77))/(JVS(179))
  XX(79) = X(79)/JVS(181)
  XX(80) = X(80)/JVS(183)
  XX(81) = X(81)/JVS(185)
  XX(82) = X(82)/JVS(187)
  XX(83) = X(83)/JVS(189)
  XX(84) = X(84)/JVS(191)
  XX(85) = X(85)/JVS(193)
  XX(86) = X(86)/JVS(195)
  XX(87) = X(87)/JVS(197)
  XX(88) = X(88)/JVS(199)
  XX(89) = X(89)/JVS(201)
  XX(90) = X(90)/JVS(203)
  XX(91) = X(91)/JVS(205)
  XX(92) = X(92)/JVS(207)
  XX(93) = X(93)/JVS(209)
  XX(94) = X(94)/JVS(211)
  XX(95) = X(95)/JVS(213)
  XX(96) = X(96)/JVS(215)
  XX(97) = X(97)/JVS(217)
  XX(98) = X(98)/JVS(219)
  XX(99) = X(99)/JVS(221)
  XX(100) = X(100)/JVS(223)
  XX(101) = X(101)/JVS(225)
  XX(102) = (X(102)-JVS(108)*XX(52)-JVS(118)*XX(56))/(JVS(227))
  XX(103) = (X(103)-JVS(114)*XX(54)-JVS(119)*XX(56))/(JVS(229))
  XX(104) = (X(104)-JVS(24)*XX(19))/(JVS(232))
  XX(105) = (X(105)-JVS(157)*XX(69)-JVS(182)*XX(79))/(JVS(238))
  XX(106) = (X(106)-JVS(63)*XX(31)-JVS(212)*XX(94))/(JVS(242))
  XX(107) = (X(107)-JVS(33)*XX(20)-JVS(196)*XX(86))/(JVS(245))
  XX(108) = (X(108)-JVS(44)*XX(23)-JVS(200)*XX(88))/(JVS(249))
  XX(109) = (X(109)-JVS(71)*XX(35)-JVS(250)*XX(108))/(JVS(254))
  XX(110) = (X(110)-JVS(57)*XX(29))/(JVS(256))
  XX(111) = (X(111)-JVS(61)*XX(30)-JVS(210)*XX(93))/(JVS(260))
  XX(112) = (X(112)-JVS(73)*XX(36)-JVS(261)*XX(111))/(JVS(265))
  XX(113) = (X(113)-JVS(65)*XX(32)-JVS(214)*XX(95))/(JVS(268))
  XX(114) = (X(114)-JVS(36)*XX(21))/(JVS(270))
  XX(115) = (X(115)-JVS(37)*XX(21)-JVS(271)*XX(114))/(JVS(274))
  XX(116) = (X(116)-JVS(105)*XX(51)-JVS(110)*XX(53)-JVS(120)*XX(56)-JVS(226)*XX(101)-JVS(228)*XX(102)-JVS(230)*XX(103))&
              &/(JVS(280))
  XX(117) = (X(117)-JVS(162)*XX(71))/(JVS(283))
  XX(118) = (X(118)-JVS(163)*XX(71)-JVS(284)*XX(117))/(JVS(288))
  XX(119) = (X(119)-JVS(164)*XX(71)-JVS(192)*XX(84)-JVS(285)*XX(117)-JVS(289)*XX(118))/(JVS(297))
  XX(120) = (X(120)-JVS(168)*XX(73)-JVS(222)*XX(99)-JVS(233)*XX(104))/(JVS(303))
  XX(121) = (X(121)-JVS(53)*XX(27)-JVS(206)*XX(91))/(JVS(307))
  XX(122) = (X(122)-JVS(55)*XX(28)-JVS(208)*XX(92))/(JVS(311))
  XX(123) = (X(123)-JVS(25)*XX(19))/(JVS(314))
  XX(124) = (X(124)-JVS(58)*XX(29)-JVS(257)*XX(110))/(JVS(320))
  XX(125) = (X(125)-JVS(38)*XX(21)-JVS(198)*XX(87)-JVS(275)*XX(115))/(JVS(327))
  XX(126) = (X(126)-JVS(99)*XX(49)-JVS(170)*XX(74))/(JVS(338))
  XX(127) = (X(127)-JVS(166)*XX(72)-JVS(220)*XX(98)-JVS(234)*XX(104)-JVS(315)*XX(123)-JVS(321)*XX(124))/(JVS(343))
  XX(128) = (X(128)-JVS(26)*XX(19))/(JVS(346))
  XX(129) = (X(129)-JVS(140)*XX(64)-JVS(328)*XX(125))/(JVS(351))
  XX(130) = (X(130)-JVS(102)*XX(50))/(JVS(356))
  XX(131) = (X(131)-JVS(40)*XX(22)-JVS(141)*XX(64)-JVS(329)*XX(125))/(JVS(361))
  XX(132) = (X(132)-JVS(41)*XX(22)-JVS(142)*XX(64)-JVS(330)*XX(125)-JVS(352)*XX(129)-JVS(362)*XX(131))/(JVS(367))
  XX(133) = (X(133)-JVS(144)*XX(65)-JVS(331)*XX(125))/(JVS(373))
  XX(134) = (X(134)-JVS(145)*XX(65)-JVS(190)*XX(83)-JVS(332)*XX(125)-JVS(374)*XX(133))/(JVS(379))
  XX(135) = (X(135)-JVS(147)*XX(66)-JVS(184)*XX(80))/(JVS(384))
  XX(136) = (X(136)-JVS(148)*XX(66)-JVS(385)*XX(135))/(JVS(391))
  XX(137) = (X(137)-JVS(137)*XX(63)-JVS(353)*XX(129)-JVS(368)*XX(132))/(JVS(397))
  XX(138) = (X(138)-JVS(50)*XX(26)-JVS(333)*XX(125))/(JVS(401))
  XX(139) = (X(139)-JVS(159)*XX(70))/(JVS(406))
  XX(140) = (X(140)-JVS(160)*XX(70)-JVS(407)*XX(139))/(JVS(411))
  XX(141) = (X(141)-JVS(138)*XX(63)-JVS(369)*XX(132)-JVS(398)*XX(137))/(JVS(419))
  XX(142) = (X(142)-JVS(51)*XX(26)-JVS(204)*XX(90)-JVS(334)*XX(125)-JVS(402)*XX(138))/(JVS(424))
  XX(143) = (X(143)-JVS(27)*XX(19)-JVS(316)*XX(123)-JVS(322)*XX(124)-JVS(347)*XX(128))/(JVS(432))
  XX(144) = (X(144)-JVS(150)*XX(67)-JVS(186)*XX(81))/(JVS(439))
  XX(145) = (X(145)-JVS(28)*XX(19)-JVS(224)*XX(100))/(JVS(447))
  XX(146) = (X(146)-JVS(46)*XX(24)-JVS(194)*XX(85))/(JVS(453))
  XX(147) = (X(147)-JVS(69)*XX(34)-JVS(111)*XX(53)-JVS(246)*XX(107)-JVS(281)*XX(116))/(JVS(463))
  XX(148) = (X(148)-JVS(29)*XX(19)-JVS(151)*XX(67)-JVS(425)*XX(142)-JVS(440)*XX(144))/(JVS(471))
  XX(149) = (X(149)-JVS(152)*XX(67)-JVS(426)*XX(142)-JVS(441)*XX(144)-JVS(472)*XX(148))/(JVS(481))
  XX(150) = (X(150)-JVS(103)*XX(50)-JVS(218)*XX(97)-JVS(357)*XX(130)-JVS(464)*XX(147))/(JVS(489))
  XX(151) = (X(151)-JVS(154)*XX(68)-JVS(358)*XX(130)-JVS(454)*XX(146)-JVS(490)*XX(150))/(JVS(495))
  XX(152) = (X(152)-JVS(100)*XX(49)-JVS(339)*XX(126)-JVS(465)*XX(147))/(JVS(506))
  XX(153) = (X(153)-JVS(251)*XX(108)-JVS(255)*XX(109)-JVS(286)*XX(117)-JVS(290)*XX(118)-JVS(298)*XX(119)-JVS(473)&
              &*XX(148)-JVS(482)*XX(149))/(JVS(529))
  XX(154) = (X(154)-JVS(59)*XX(29)-JVS(112)*XX(53)-JVS(135)*XX(62)-JVS(175)*XX(76)-JVS(180)*XX(78)-JVS(231)*XX(103)&
              &-JVS(282)*XX(116)-JVS(317)*XX(123)-JVS(323)*XX(124)-JVS(344)*XX(127)-JVS(433)*XX(143)-JVS(466)*XX(147)&
              &-JVS(530)*XX(153))/(JVS(556))
  XX(155) = (X(155)-JVS(155)*XX(68)-JVS(188)*XX(82)-JVS(348)*XX(128)-JVS(434)*XX(143)-JVS(455)*XX(146)-JVS(491)*XX(150)&
              &-JVS(496)*XX(151)-JVS(507)*XX(152)-JVS(531)*XX(153)-JVS(557)*XX(154))/(JVS(585))
  XX(156) = (X(156)-JVS(87)*XX(43)-JVS(247)*XX(107)-JVS(258)*XX(110)-JVS(262)*XX(111)-JVS(266)*XX(112)-JVS(272)*XX(114)&
              &-JVS(276)*XX(115)-JVS(308)*XX(121)-JVS(312)*XX(122)-JVS(324)*XX(124)-JVS(335)*XX(125)-JVS(340)*XX(126)&
              &-JVS(349)*XX(128)-JVS(354)*XX(129)-JVS(359)*XX(130)-JVS(363)*XX(131)-JVS(370)*XX(132)-JVS(375)*XX(133)&
              &-JVS(380)*XX(134)-JVS(386)*XX(135)-JVS(392)*XX(136)-JVS(399)*XX(137)-JVS(403)*XX(138)-JVS(408)*XX(139)&
              &-JVS(412)*XX(140)-JVS(420)*XX(141)-JVS(427)*XX(142)-JVS(448)*XX(145)-JVS(456)*XX(146)-JVS(467)*XX(147)&
              &-JVS(483)*XX(149)-JVS(492)*XX(150)-JVS(497)*XX(151)-JVS(508)*XX(152)-JVS(532)*XX(153)-JVS(558)*XX(154)&
              &-JVS(586)*XX(155))/(JVS(633))
  XX(157) = (X(157)-JVS(30)*XX(19)-JVS(235)*XX(104)-JVS(291)*XX(118)-JVS(299)*XX(119)-JVS(304)*XX(120)-JVS(309)*XX(121)&
              &-JVS(313)*XX(122)-JVS(345)*XX(127)-JVS(350)*XX(128)-JVS(435)*XX(143)-JVS(442)*XX(144)-JVS(449)*XX(145)&
              &-JVS(457)*XX(146)-JVS(468)*XX(147)-JVS(474)*XX(148)-JVS(484)*XX(149)-JVS(493)*XX(150)-JVS(509)*XX(152)&
              &-JVS(533)*XX(153)-JVS(559)*XX(154)-JVS(587)*XX(155)-JVS(634)*XX(156))/(JVS(658))
  XX(158) = (X(158)-JVS(48)*XX(25)-JVS(202)*XX(89)-JVS(450)*XX(145)-JVS(469)*XX(147)-JVS(498)*XX(151)-JVS(510)*XX(152)&
              &-JVS(534)*XX(153)-JVS(560)*XX(154)-JVS(588)*XX(155)-JVS(635)*XX(156)-JVS(659)*XX(157))/(JVS(670))
  XX(159) = (X(159)-JVS(97)*XX(48)-JVS(216)*XX(96)-JVS(360)*XX(130)-JVS(451)*XX(145)-JVS(494)*XX(150)-JVS(499)*XX(151)&
              &-JVS(511)*XX(152)-JVS(535)*XX(153)-JVS(561)*XX(154)-JVS(589)*XX(155)-JVS(636)*XX(156)-JVS(660)*XX(157)&
              &-JVS(671)*XX(158))/(JVS(686))
  XX(160) = (X(160)-JVS(31)*XX(19)-JVS(292)*XX(118)-JVS(300)*XX(119)-JVS(443)*XX(144)-JVS(475)*XX(148)-JVS(485)*XX(149)&
              &-JVS(512)*XX(152)-JVS(536)*XX(153)-JVS(562)*XX(154)-JVS(590)*XX(155)-JVS(637)*XX(156)-JVS(661)*XX(157)&
              &-JVS(672)*XX(158)-JVS(687)*XX(159))/(JVS(701))
  XX(161) = (X(161)-JVS(42)*XX(22)-JVS(67)*XX(33)-JVS(236)*XX(104)-JVS(239)*XX(105)-JVS(243)*XX(106)-JVS(269)*XX(113)&
              &-JVS(287)*XX(117)-JVS(293)*XX(118)-JVS(301)*XX(119)-JVS(336)*XX(125)-JVS(341)*XX(126)-JVS(355)*XX(129)&
              &-JVS(364)*XX(131)-JVS(371)*XX(132)-JVS(376)*XX(133)-JVS(381)*XX(134)-JVS(387)*XX(135)-JVS(393)*XX(136)&
              &-JVS(400)*XX(137)-JVS(404)*XX(138)-JVS(409)*XX(139)-JVS(413)*XX(140)-JVS(421)*XX(141)-JVS(428)*XX(142)&
              &-JVS(436)*XX(143)-JVS(458)*XX(146)-JVS(476)*XX(148)-JVS(486)*XX(149)-JVS(500)*XX(151)-JVS(513)*XX(152)&
              &-JVS(537)*XX(153)-JVS(563)*XX(154)-JVS(591)*XX(155)-JVS(638)*XX(156)-JVS(662)*XX(157)-JVS(673)*XX(158)&
              &-JVS(688)*XX(159)-JVS(702)*XX(160))/(JVS(738))
  XX(161) = XX(161)
  XX(160) = XX(160)-JVS(737)*XX(161)
  XX(159) = XX(159)-JVS(700)*XX(160)-JVS(736)*XX(161)
  XX(158) = XX(158)-JVS(685)*XX(159)-JVS(699)*XX(160)-JVS(735)*XX(161)
  XX(157) = XX(157)-JVS(669)*XX(158)-JVS(684)*XX(159)-JVS(698)*XX(160)-JVS(734)*XX(161)
  XX(156) = XX(156)-JVS(657)*XX(157)-JVS(668)*XX(158)-JVS(683)*XX(159)-JVS(697)*XX(160)-JVS(733)*XX(161)
  XX(155) = XX(155)-JVS(632)*XX(156)-JVS(656)*XX(157)-JVS(667)*XX(158)-JVS(682)*XX(159)-JVS(696)*XX(160)-JVS(732)&
              &*XX(161)
  XX(154) = XX(154)-JVS(584)*XX(155)-JVS(631)*XX(156)-JVS(655)*XX(157)-JVS(681)*XX(159)-JVS(695)*XX(160)-JVS(731)&
              &*XX(161)
  XX(153) = XX(153)-JVS(583)*XX(155)-JVS(630)*XX(156)-JVS(654)*XX(157)-JVS(694)*XX(160)-JVS(730)*XX(161)
  XX(152) = XX(152)-JVS(528)*XX(153)-JVS(555)*XX(154)-JVS(629)*XX(156)-JVS(666)*XX(158)-JVS(680)*XX(159)-JVS(729)&
              &*XX(161)
  XX(151) = XX(151)-JVS(505)*XX(152)-JVS(527)*XX(153)-JVS(554)*XX(154)-JVS(582)*XX(155)-JVS(628)*XX(156)-JVS(653)&
              &*XX(157)-JVS(665)*XX(158)-JVS(679)*XX(159)-JVS(728)*XX(161)
  XX(150) = XX(150)-JVS(504)*XX(152)-JVS(526)*XX(153)-JVS(553)*XX(154)-JVS(581)*XX(155)-JVS(627)*XX(156)-JVS(652)&
              &*XX(157)-JVS(678)*XX(159)-JVS(727)*XX(161)
  XX(149) = XX(149)-JVS(525)*XX(153)-JVS(580)*XX(155)-JVS(626)*XX(156)-JVS(651)*XX(157)-JVS(693)*XX(160)-JVS(726)&
              &*XX(161)
  XX(148) = XX(148)-JVS(480)*XX(149)-JVS(524)*XX(153)-JVS(579)*XX(155)-JVS(625)*XX(156)-JVS(650)*XX(157)-JVS(692)&
              &*XX(160)-JVS(725)*XX(161)
  XX(147) = XX(147)-JVS(523)*XX(153)-JVS(552)*XX(154)-JVS(624)*XX(156)-JVS(677)*XX(159)-JVS(724)*XX(161)
  XX(146) = XX(146)-JVS(551)*XX(154)-JVS(578)*XX(155)-JVS(623)*XX(156)-JVS(649)*XX(157)-JVS(723)*XX(161)
  XX(145) = XX(145)-JVS(503)*XX(152)-JVS(550)*XX(154)-JVS(622)*XX(156)-JVS(648)*XX(157)-JVS(664)*XX(158)-JVS(676)&
              &*XX(159)
  XX(144) = XX(144)-JVS(470)*XX(148)-JVS(479)*XX(149)-JVS(522)*XX(153)-JVS(647)*XX(157)-JVS(691)*XX(160)
  XX(143) = XX(143)-JVS(549)*XX(154)-JVS(577)*XX(155)-JVS(621)*XX(156)-JVS(646)*XX(157)-JVS(722)*XX(161)
  XX(142) = XX(142)-JVS(521)*XX(153)-JVS(576)*XX(155)-JVS(620)*XX(156)-JVS(721)*XX(161)
  XX(141) = XX(141)-JVS(520)*XX(153)-JVS(575)*XX(155)-JVS(619)*XX(156)-JVS(720)*XX(161)
  XX(140) = XX(140)-JVS(418)*XX(141)-JVS(574)*XX(155)-JVS(618)*XX(156)-JVS(719)*XX(161)
  XX(139) = XX(139)-JVS(410)*XX(140)-JVS(417)*XX(141)-JVS(573)*XX(155)-JVS(617)*XX(156)-JVS(718)*XX(161)
  XX(138) = XX(138)-JVS(423)*XX(142)-JVS(519)*XX(153)-JVS(572)*XX(155)-JVS(616)*XX(156)-JVS(717)*XX(161)
  XX(137) = XX(137)-JVS(416)*XX(141)-JVS(518)*XX(153)-JVS(571)*XX(155)-JVS(615)*XX(156)-JVS(716)*XX(161)
  XX(136) = XX(136)-JVS(396)*XX(137)-JVS(415)*XX(141)-JVS(478)*XX(149)-JVS(614)*XX(156)-JVS(715)*XX(161)
  XX(135) = XX(135)-JVS(390)*XX(136)-JVS(395)*XX(137)-JVS(477)*XX(149)-JVS(613)*XX(156)-JVS(714)*XX(161)
  XX(134) = XX(134)-JVS(383)*XX(135)-JVS(389)*XX(136)-JVS(570)*XX(155)-JVS(612)*XX(156)-JVS(713)*XX(161)
  XX(133) = XX(133)-JVS(378)*XX(134)-JVS(388)*XX(136)-JVS(569)*XX(155)-JVS(611)*XX(156)-JVS(712)*XX(161)
  XX(132) = XX(132)-JVS(517)*XX(153)-JVS(610)*XX(156)-JVS(711)*XX(161)
  XX(131) = XX(131)-JVS(366)*XX(132)-JVS(516)*XX(153)-JVS(609)*XX(156)-JVS(710)*XX(161)
  XX(130) = XX(130)-JVS(488)*XX(150)-JVS(502)*XX(152)-JVS(608)*XX(156)
  XX(129) = XX(129)-JVS(365)*XX(132)-JVS(607)*XX(156)-JVS(709)*XX(161)
  XX(128) = XX(128)-JVS(548)*XX(154)-JVS(606)*XX(156)-JVS(645)*XX(157)
  XX(127) = XX(127)-JVS(431)*XX(143)-JVS(462)*XX(147)-JVS(547)*XX(154)-JVS(605)*XX(156)-JVS(644)*XX(157)-JVS(708)&
              &*XX(161)
  XX(126) = XX(126)-JVS(501)*XX(152)-JVS(604)*XX(156)-JVS(675)*XX(159)-JVS(707)*XX(161)
  XX(125) = XX(125)-JVS(603)*XX(156)
  XX(124) = XX(124)-JVS(546)*XX(154)-JVS(602)*XX(156)
  XX(123) = XX(123)-JVS(319)*XX(124)-JVS(706)*XX(161)
  XX(122) = XX(122)-JVS(545)*XX(154)-JVS(601)*XX(156)-JVS(643)*XX(157)
  XX(121) = XX(121)-JVS(544)*XX(154)-JVS(600)*XX(156)-JVS(642)*XX(157)
  XX(120) = XX(120)-JVS(306)*XX(121)-JVS(430)*XX(143)-JVS(438)*XX(144)-JVS(568)*XX(155)-JVS(641)*XX(157)
  XX(119) = XX(119)
  XX(118) = XX(118)-JVS(296)*XX(119)
  XX(117) = XX(117)-JVS(295)*XX(119)
  XX(116) = XX(116)-JVS(461)*XX(147)-JVS(543)*XX(154)
  XX(115) = XX(115)-JVS(326)*XX(125)-JVS(599)*XX(156)
  XX(114) = XX(114)-JVS(273)*XX(115)-JVS(598)*XX(156)
  XX(113) = XX(113)-JVS(372)*XX(133)-JVS(414)*XX(141)-JVS(567)*XX(155)-JVS(705)*XX(161)
  XX(112) = XX(112)-JVS(597)*XX(156)
  XX(111) = XX(111)-JVS(264)*XX(112)-JVS(596)*XX(156)
  XX(110) = XX(110)-JVS(318)*XX(124)-JVS(595)*XX(156)
  XX(109) = XX(109)-JVS(515)*XX(153)
  XX(108) = XX(108)-JVS(253)*XX(109)-JVS(514)*XX(153)
  XX(107) = XX(107)-JVS(460)*XX(147)-JVS(594)*XX(156)
  XX(106) = XX(106)-JVS(394)*XX(137)-JVS(566)*XX(155)-JVS(704)*XX(161)
  XX(105) = XX(105)-JVS(241)*XX(106)-JVS(405)*XX(139)-JVS(565)*XX(155)-JVS(703)*XX(161)
  XX(104) = XX(104)-JVS(429)*XX(143)
  XX(103) = XX(103)-JVS(279)*XX(116)
  XX(102) = XX(102)-JVS(278)*XX(116)
  XX(101) = XX(101)-JVS(277)*XX(116)
  XX(100) = XX(100)-JVS(446)*XX(145)
  XX(99) = XX(99)-JVS(302)*XX(120)
  XX(98) = XX(98)-JVS(342)*XX(127)
  XX(97) = XX(97)-JVS(487)*XX(150)
  XX(96) = XX(96)-JVS(674)*XX(159)
  XX(95) = XX(95)-JVS(267)*XX(113)
  XX(94) = XX(94)-JVS(240)*XX(106)
  XX(93) = XX(93)-JVS(259)*XX(111)
  XX(92) = XX(92)-JVS(310)*XX(122)
  XX(91) = XX(91)-JVS(305)*XX(121)
  XX(90) = XX(90)-JVS(422)*XX(142)
  XX(89) = XX(89)-JVS(663)*XX(158)
  XX(88) = XX(88)-JVS(248)*XX(108)
  XX(87) = XX(87)-JVS(325)*XX(125)
  XX(86) = XX(86)-JVS(244)*XX(107)
  XX(85) = XX(85)-JVS(452)*XX(146)
  XX(84) = XX(84)-JVS(294)*XX(119)
  XX(83) = XX(83)-JVS(377)*XX(134)
  XX(82) = XX(82)-JVS(564)*XX(155)
  XX(81) = XX(81)-JVS(437)*XX(144)
  XX(80) = XX(80)-JVS(382)*XX(135)
  XX(79) = XX(79)-JVS(237)*XX(105)
  XX(78) = XX(78)-JVS(542)*XX(154)
  XX(77) = XX(77)-JVS(178)*XX(78)
  XX(76) = XX(76)-JVS(541)*XX(154)
  XX(75) = XX(75)-JVS(173)*XX(76)
  XX(74) = XX(74)-JVS(337)*XX(126)
  XX(73) = XX(73)
  XX(72) = XX(72)
  XX(71) = XX(71)
  XX(70) = XX(70)
  XX(69) = XX(69)
  XX(68) = XX(68)
  XX(67) = XX(67)
  XX(66) = XX(66)
  XX(65) = XX(65)
  XX(64) = XX(64)
  XX(63) = XX(63)
  XX(62) = XX(62)
  XX(61) = XX(61)
  XX(60) = XX(60)
  XX(59) = XX(59)
  XX(58) = XX(58)
  XX(57) = XX(57)
  XX(56) = XX(56)
  XX(55) = XX(55)
  XX(54) = XX(54)
  XX(53) = XX(53)
  XX(52) = XX(52)
  XX(51) = XX(51)
  XX(50) = XX(50)
  XX(49) = XX(49)
  XX(48) = XX(48)
  XX(47) = XX(47)
  XX(46) = XX(46)
  XX(45) = XX(45)
  XX(44) = XX(44)
  XX(43) = XX(43)
  XX(42) = XX(42)
  XX(41) = XX(41)
  XX(40) = XX(40)
  XX(39) = XX(39)
  XX(38) = XX(38)
  XX(37) = XX(37)
  XX(36) = XX(36)
  XX(35) = XX(35)
  XX(34) = XX(34)
  XX(33) = XX(33)
  XX(32) = XX(32)
  XX(31) = XX(31)
  XX(30) = XX(30)
  XX(29) = XX(29)
  XX(28) = XX(28)
  XX(27) = XX(27)
  XX(26) = XX(26)
  XX(25) = XX(25)
  XX(24) = XX(24)
  XX(23) = XX(23)
  XX(22) = XX(22)
  XX(21) = XX(21)
  XX(20) = XX(20)
  XX(19) = XX(19)
  XX(18) = XX(18)-JVS(94)*XX(47)
  XX(17) = XX(17)-JVS(92)*XX(46)
  XX(16) = XX(16)-JVS(90)*XX(45)
  XX(15) = XX(15)-JVS(88)*XX(44)
  XX(14) = XX(14)-JVS(78)*XX(39)-JVS(445)*XX(145)-JVS(540)*XX(154)-JVS(593)*XX(156)-JVS(640)*XX(157)-JVS(690)*XX(160)
  XX(13) = XX(13)-JVS(76)*XX(38)-JVS(444)*XX(145)-JVS(539)*XX(154)-JVS(592)*XX(156)-JVS(639)*XX(157)-JVS(689)*XX(160)
  XX(12) = XX(12)-JVS(84)*XX(42)
  XX(11) = XX(11)-JVS(82)*XX(41)
  XX(10) = XX(10)-JVS(80)*XX(40)
  XX(9) = XX(9)-JVS(74)*XX(37)
  XX(8) = XX(8)-JVS(9)*XX(9)
  XX(7) = XX(7)-JVS(14)*XX(12)
  XX(6) = XX(6)-JVS(11)*XX(10)
  XX(5) = XX(5)-JVS(538)*XX(154)
  XX(4) = XX(4)-JVS(21)*XX(18)
  XX(3) = XX(3)-JVS(263)*XX(112)
  XX(2) = XX(2)-JVS(252)*XX(109)
  XX(1) = XX(1)-JVS(459)*XX(147)
      
END SUBROUTINE KppSolveTR

! End of KppSolveTR function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! BLAS_UTIL - BLAS-LIKE utility functions
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!--------------------------------------------------------------
!
! BLAS/LAPACK-like subroutines used by the integration algorithms
! It is recommended to replace them by calls to the optimized
!      BLAS/LAPACK library for your machine
!
!  (C) Adrian Sandu, Aug. 2004
!      Virginia Polytechnic Institute and State University
!--------------------------------------------------------------


!--------------------------------------------------------------
      SUBROUTINE WCOPY(N,X,incX,Y,incY)
!--------------------------------------------------------------
!     copies a vector, x, to a vector, y:  y <- x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
!--------------------------------------------------------------
!     USE aqchem_Precision
      
      INTEGER  :: i,incX,incY,M,MP1,N
      REAL(kind=dp) :: X(N),Y(N)

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = X(i)
        END DO
        IF( N .LT. 8 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,8
        Y(i) = X(i)
        Y(i + 1) = X(i + 1)
        Y(i + 2) = X(i + 2)
        Y(i + 3) = X(i + 3)
        Y(i + 4) = X(i + 4)
        Y(i + 5) = X(i + 5)
        Y(i + 6) = X(i + 6)
        Y(i + 7) = X(i + 7)
      END DO

      END SUBROUTINE WCOPY


!--------------------------------------------------------------
      SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
!--------------------------------------------------------------
!     constant times a vector plus a vector: y <- y + Alpha*x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
!--------------------------------------------------------------

      INTEGER  :: i,incX,incY,M,MP1,N
      REAL(kind=dp) :: X(N),Y(N),Alpha
      REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp

      IF (Alpha .EQ. ZERO) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,4)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = Y(i) + Alpha*X(i)
        END DO
        IF( N .LT. 4 ) RETURN
      END IF
      MP1 = M + 1
      DO i = MP1,N,4
        Y(i) = Y(i) + Alpha*X(i)
        Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
        Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
        Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
      END DO
      
      END SUBROUTINE WAXPY



!--------------------------------------------------------------
      SUBROUTINE WSCAL(N,Alpha,X,incX)
!--------------------------------------------------------------
!     constant times a vector: x(1:N) <- Alpha*x(1:N) 
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
!--------------------------------------------------------------

      INTEGER  :: i,incX,M,MP1,N
      REAL(kind=dp)  :: X(N),Alpha
      REAL(kind=dp), PARAMETER  :: ZERO=0.0_dp, ONE=1.0_dp

      IF (Alpha .EQ. ONE) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,5)
      IF( M .NE. 0 ) THEN
        IF (Alpha .EQ. (-ONE)) THEN
          DO i = 1,M
            X(i) = -X(i)
          END DO
        ELSEIF (Alpha .EQ. ZERO) THEN
          DO i = 1,M
            X(i) = ZERO
          END DO
        ELSE
          DO i = 1,M
            X(i) = Alpha*X(i)
          END DO
        END IF
        IF( N .LT. 5 ) RETURN
      END IF
      MP1 = M + 1
      IF (Alpha .EQ. (-ONE)) THEN
        DO i = MP1,N,5
          X(i)     = -X(i)
          X(i + 1) = -X(i + 1)
          X(i + 2) = -X(i + 2)
          X(i + 3) = -X(i + 3)
          X(i + 4) = -X(i + 4)
        END DO
      ELSEIF (Alpha .EQ. ZERO) THEN
        DO i = MP1,N,5
          X(i)     = ZERO
          X(i + 1) = ZERO
          X(i + 2) = ZERO
          X(i + 3) = ZERO
          X(i + 4) = ZERO
        END DO
      ELSE
        DO i = MP1,N,5
          X(i)     = Alpha*X(i)
          X(i + 1) = Alpha*X(i + 1)
          X(i + 2) = Alpha*X(i + 2)
          X(i + 3) = Alpha*X(i + 3)
          X(i + 4) = Alpha*X(i + 4)
        END DO
      END IF

      END SUBROUTINE WSCAL

!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION WLAMCH( C )
!--------------------------------------------------------------
!     returns epsilon machine
!     after LAPACK
!     replace this by the function from the optimized LAPACK implementation:
!          CALL SLAMCH('E') or CALL DLAMCH('E')
!--------------------------------------------------------------
!      USE aqchem_Precision

      CHARACTER ::  C
      INTEGER    :: i
      REAL(kind=dp), SAVE  ::  Eps
      REAL(kind=dp)  ::  Suma
      REAL(kind=dp), PARAMETER  ::  ONE=1.0_dp, HALF=0.5_dp
      LOGICAL, SAVE   ::  First=.TRUE.
      
      IF (First) THEN
        First = .FALSE.
        Eps = HALF**(16)
        DO i = 17, 80
          Eps = Eps*HALF
          CALL WLAMCH_ADD(ONE,Eps,Suma)
          IF (Suma.LE.ONE) GOTO 10
        END DO
        PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
        RETURN
10      Eps = Eps*2
        i = i-1      
      END IF

      WLAMCH = Eps

      END FUNCTION WLAMCH
     
      SUBROUTINE WLAMCH_ADD( A, B, Suma )
!      USE aqchem_Precision
      
      REAL(kind=dp) A, B, Suma
      Suma = A + B

      END SUBROUTINE WLAMCH_ADD
!--------------------------------------------------------------


!--------------------------------------------------------------
      SUBROUTINE SET2ZERO(N,Y)
!--------------------------------------------------------------
!     copies zeros into the vector y:  y <- 0
!     after BLAS
!--------------------------------------------------------------
      
      INTEGER ::  i,M,MP1,N
      REAL(kind=dp) ::  Y(N)
      REAL(kind=dp), PARAMETER :: ZERO = 0.0d0

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = ZERO
        END DO
        IF( N .LT. 8 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,8
        Y(i)     = ZERO
        Y(i + 1) = ZERO
        Y(i + 2) = ZERO
        Y(i + 3) = ZERO
        Y(i + 4) = ZERO
        Y(i + 5) = ZERO
        Y(i + 6) = ZERO
        Y(i + 7) = ZERO
      END DO

      END SUBROUTINE SET2ZERO


!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION WDOT (N, DX, incX, DY, incY) 
!--------------------------------------------------------------
!     dot produce: wdot = x(1:N)*y(1:N) 
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
!--------------------------------------------------------------
!      USE messy_mecca_kpp_Precision
!--------------------------------------------------------------
      IMPLICIT NONE
      INTEGER :: N, incX, incY
      REAL(kind=dp) :: DX(N), DY(N) 

      INTEGER :: i, IX, IY, M, MP1, NS
                                 
      WDOT = 0.0D0 
      IF (N .LE. 0) RETURN 
      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
!                                                                       
!     Code for unequal or nonpositive increments.                       
!                                                                       
    5 IX = 1 
      IY = 1 
      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
      DO i = 1,N 
        WDOT = WDOT + DX(IX)*DY(IY) 
        IX = IX + incX 
        IY = IY + incY 
      END DO 
      RETURN 
!                                                                       
!     Code for both increments equal to 1.                              
!                                                                       
!     Clean-up loop so remaining vector length is a multiple of 5.      
!                                                                       
   20 M = MOD(N,5) 
      IF (M .EQ. 0) GO TO 40 
      DO i = 1,M 
         WDOT = WDOT + DX(i)*DY(i) 
      END DO 
      IF (N .LT. 5) RETURN 
   40 MP1 = M + 1 
      DO i = MP1,N,5 
          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
      END DO 
      RETURN 
!                                                                       
!     Code for equal, positive, non-unit increments.                    
!                                                                       
   60 NS = N*incX 
      DO i = 1,NS,incX 
        WDOT = WDOT + DX(i)*DY(i) 
      END DO 

      END FUNCTION WDOT                                          


!--------------------------------------------------------------
      SUBROUTINE WADD(N,X,Y,Z)
!--------------------------------------------------------------
!     adds two vectors: z <- x + y
!     BLAS - like
!--------------------------------------------------------------
!     USE aqchem_Precision
      
      INTEGER :: i, M, MP1, N
      REAL(kind=dp) :: X(N),Y(N),Z(N)

      IF (N.LE.0) RETURN

      M = MOD(N,5)
      IF( M /= 0 ) THEN
         DO i = 1,M
            Z(i) = X(i) + Y(i)
         END DO
         IF( N < 5 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,5
         Z(i)     = X(i)     + Y(i)
         Z(i + 1) = X(i + 1) + Y(i + 1)
         Z(i + 2) = X(i + 2) + Y(i + 2)
         Z(i + 3) = X(i + 3) + Y(i + 3)
         Z(i + 4) = X(i + 4) + Y(i + 4)
      END DO

      END SUBROUTINE WADD
      
      
      
!--------------------------------------------------------------
      SUBROUTINE WGEFA(N,A,Ipvt,info)
!--------------------------------------------------------------
!     WGEFA FACTORS THE MATRIX A (N,N) BY
!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
!     LINPACK - LIKE 
!--------------------------------------------------------------
!
      INTEGER       :: N,Ipvt(N),info
      REAL(kind=dp) :: A(N,N)
      REAL(kind=dp) :: t, dmax, da
      INTEGER       :: j,k,l
      REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0

      info = 0

size: IF (n > 1) THEN
      
col:  DO k = 1, n-1

!        find l = pivot index
!        l = idamax(n-k+1,A(k,k),1) + k - 1
         l = k; dmax = abs(A(k,k))
         DO j = k+1,n
            da = ABS(A(j,k))
            IF (da > dmax) THEN
              l = j; dmax = da
            END IF
         END DO
         Ipvt(k) = l

!        zero pivot implies this column already triangularized
         IF (ABS(A(l,k)) < TINY(ZERO)) THEN
            info = k
            return
         ELSE   
            IF (l /= k) THEN
               t = A(l,k); A(l,k) = A(k,k); A(k,k) = t
            END IF
            t = -ONE/A(k,k)
            CALL WSCAL(n-k,t,A(k+1,k),1)
            DO j = k+1, n
               t = A(l,j)
               IF (l /= k) THEN
                  A(l,j) = A(k,j); A(k,j) = t
               END IF
               CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1)
            END DO         
         END IF
         
       END DO col
       
      END IF size
      
      Ipvt(N) = N
      IF (ABS(A(N,N)) == ZERO) info = N
      
      END SUBROUTINE WGEFA


!--------------------------------------------------------------
      SUBROUTINE WGESL(Trans,N,A,Ipvt,b)
!--------------------------------------------------------------
!     WGESL solves the system
!     a * x = b  or  trans(a) * x = b
!     using the factors computed by WGEFA.
!
!     Trans      = 'N'   to solve  A*x = b ,
!                = 'T'   to solve  transpose(A)*x = b
!     LINPACK - LIKE 
!--------------------------------------------------------------

      INTEGER       :: N,Ipvt(N)
      CHARACTER     :: trans
      REAL(kind=dp) :: A(N,N),b(N)
      REAL(kind=dp) :: t
      INTEGER       :: k,kb,l

      
      SELECT CASE (Trans)

      CASE ('n','N')  !  Solve  A * x = b

!        first solve  L*y = b
         IF (n >= 2) THEN
          DO k = 1, n-1
            l = Ipvt(k)
            t = b(l)
            IF (l /= k) THEN
               b(l) = b(k)
               b(k) = t
            END IF
            CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1)
          END DO
         END IF
!        now solve  U*x = y
         DO kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            CALL WAXPY(k-1,t,a(1,k),1,b(1),1)
         END DO
      
      CASE ('t','T')  !  Solve transpose(A) * x = b

!        first solve  trans(U)*y = b
         DO k = 1, n
            t = WDOT(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
         END DO
!        now solve trans(L)*x = y
         IF (n >= 2) THEN
         DO kb = 1, n-1
            k = n - kb
            b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
            l = Ipvt(k)
            IF (l /= k) THEN
               t = b(l); b(l) = b(k); b(k) = t
            END IF
         END DO
         END IF
   
      END SELECT

      END SUBROUTINE WGESL
! End of BLAS_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



END MODULE aqchem_LinearAlgebra

